function dL = derivativeL(NUM_RATE,Time,jumpR1, jumpR2, q, cond)

% Time = temp(k,1);
% jumpR1 = temp(k,2)
% jumpR2 = temp(k,3)

dL = sym('A%d%d',[8,8]);
dL(:,:) = 0;

switch cond
    
    case 'normal'
        for para_i = 1: NUM_RATE
            for para_j = 1:NUM_RATE+1
                if para_i ~= para_j
                    dP = derivativeP(para_i, para_j,Time,q);
                    p = expm(q*Time);
                    dL(para_i,para_j) = 1/p(jumpR1,jumpR2)* dP(jumpR1,jumpR2);
                end
            end
        end
        dL = double(dL);
        
    case 'cencor'
        for para_i = 1: NUM_RATE
            for para_j = 1:NUM_RATE+1
                if para_i ~= para_j
                    
                    dP = derivativeP(para_i, para_j,Time,q);  
                    
                    Pij = 0;
                    for l = 1:NUM_RATE
                        p = expm(q*Time);
                        Pij = Pij + p(jumpR2,l);                        
                    end
                    
                    dP_sum = 0;
                    for l = 1:NUM_RATE
                        dP_sum = dP_sum+ dP(jumpR2,l);
                    end
                    
                    dL(para_i,para_j) = 1/Pij* dP_sum;
                end
            end
        end
        dL = double(dL);
        
    case 'default'
        for para_i = 1: NUM_RATE
            for para_j = 1:NUM_RATE+1
                if para_i ~= para_j
                    
                    dP = derivativeP(para_i, para_j,Time,q);                    
                    Pij = 0;
                    t1 = (Time - MIN_TIME_UNIT);
                    t2 =  MIN_TIME_UNIT;
                    p1 = expm(q*t1);
                    p2 = expm(q*t2);
                    for l = 1:NUM_RATE                       
                        a = p1(jumpR1,l);
                        b = p2(l,jumpR2);
                        Pij = Pij + a*b;
                    end
                    
                    dP_sum = 0;
                    for l = 1:NUM_RATE
                        dP_sum = dP_sum + dP(jumpR1,l)* p2(l,NUM_RATE+1) * t1* t2;                        
                    end
                    
                    dL(para_i,para_j) = 1/Pij* dP_sum;
                end
            end
        end
        dL = double(dL);
        
end




end




% Derivatives of P respect to parameters
function dP = derivativeP(para_i, para_j,Time,q)

[A,D] = eig(q);

dQ = derivativeQ(para_i, para_j);

% kxk matrix G = inv(A) * dQ *A
G = double(A\dQ* A);  % convert from symbolic to double

V = zeros(8,8);
for i = 1:7
    for j = 1:7
        
        if i == j
            V(i,j) = G(i,j)* Time * exp(D(i,j)*Time);
        else
            V(i,j) = G(i,j)* (exp(D(i,i)*Time) - exp(D(j,j)*Time))/(D(i,i) - D(j,j));
        end
                
    end
end

% kxk matrix dP = A * V * inv(A)
dP = real(A*V/A);

end

% Derivatives of Q respect to parameters
function dQ = derivativeQ(para_i, para_j)
Q = sym('A%d%d', [8 8]);
for i = 1:8
    Q(i,i) = 0;
end
Q(8,:) = 0;

for i = 1: 7
    for j = 1: 8
        if j ~= i
            Q(i,i) = Q(i,i)-Q(i,j);
        end
    end
end

dQ = diff(Q,Q(para_i, para_j));
end